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ABSTRACT; 

This  paper  will  address  analytical  research  work  performed  on  impact  location  in  an 
anisotropic  multilayered  thermal  protection  structure.  Hie  propagation  of  elastodynamic 
waves  through  the  multilayer  structure  and  the  possibility  that  the  layers  are  transversely 
isotropic  are  two  complications  that  arc  addressed.  The  method  used  is  an  extension  of 
classical  triangulation  based  on  the  use  of  the  fixed  geometry  of  die  multilayer  plate 
structure,  qua  si -longitudinal  modes,  and  the  generalized  Snell’s  law  for  transversely 
isotropic  materials.  The  impact  localization  problem  is  recast  as  a  minimization  problem 
whose  objective  function  is  the  sum  of  the  distances  between  all  combination  pairs  of 
constant  time  difference  curves  for  the  sensor  pairs.  An  initial  estimate  fs  made  for  the 
parameter  values  for  the  collection  of  constant  time  difference  curves.  The 
Newton-Kantorovich  algorithm  is  then  used  to  generate  a  sequence  of  iterations  which 
converge  to  the  minimum  of  the  objective  function,  This  generates  a  point  on  each  of  the 
collection  of  curves.  The  centroid  of  this  point  set  is  calculated  and  used  for  the 
approximation  for  the  impact  point. 


INTRODUCTION 

Thermal  protection  systems  on  aerospace  vehicles  are  exposed  to  possible  damage  by 
impact  from  debris,  maintenance  equipment,  etc.  The  ability  would  exist  in  an  ideal  system 
to  detect  that  an  impact  event  has  occurred,  where  the  impact  has  occurred,  quantify  the 
impact  energy,  and  determine  the  extent  to  which  damage  has  occurred,  A  passive  system  to 
accomplish  these  goals  using  the  acoustic  energy  generated  by  the  impact  itself  would  be 
desirable.  The  subject  of  the  current  work  is  the  second  level:  impact  localization. 
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A  thermal  protection  system  (TPS)  is  usually  composed  of  an  outer  refractory  layer  with 
intermediate  bonding  or  strain  isolation  layers  which  tie  it  to  the  vehicle  outer  surface 


Figure  L  Multilayer  Transversely  Isotropic  Flatc(side).  Figure  2.  Multilayer  Transversely  Isotropic  Pia(c(bottom). 

Due  to  the  extreme  environment  on  the  exterior  of  the  TPS  the  use  of  backside  sensing  is 
considered  to  be  practical.  Thus  the  signals  that  will  be  detected  from  an  impact  will  have 
to  pass  thru  this  layered  structure.  Some  of  the  materials  in  this  structure  may  also  be 
anisotropic.  The  ceramic  foam  tile  used  in  some  TPS  is  transversely  isotropic.  This  leads 
to  the  problem  illustrated  in  Figure  1  and  2. 

It  is  assumed  that  the  impact  generates  signals  that  propagate  out  from  the  location  of 
impact  through  the  outermost  layer  media.  The  impact  signals  detected  by  the  sensors  have  a 
propagation  delay  and  the  actual  time  an  impact  has  occurred  is  unknown.  This  problem  is 
similar  in  character  to  the  classic  location  problem  in  an  isotropic  medium  where  two 
propagation  observations  are  required  if  the  time  of  signal  initiation  is  know'  or  three  are 
required  if  the  initiation  time  of  the  signal  is  unknown.  In  this  case,  the  propagation  time  is 
directly  proportional  to  the  propagation  distance  and  may  be  direc  tly  solved  for.  The  layered 
structure  and  anisotropic  material  properties  make  the  explicit  analytic  statement  of  the 
relationship  between  propagation  time  and  distance  difficult,  if  not  impossible.  This  leads  to 
a  numerical  approach  for  finding  a  solution. 

This  situation  may  be  viewed  as  a  Model  Based  Inverse  problem.  If  the  location  of 
impact  and  initial  propagation  direction  is  specified  then  the  time  to  propagate  through  the 
layer  assembly  and  exit  point  on  the  backside  surface  may  be  calculated.  This  would  be  the 
forward  problem.  The  inverse  problem  consists  of  knowledge  of  the  relative  arrival  times  at 
various  monitoring  points  on  the  backside  surface  and  from  this  information  finding  the 
impact  coordinates. 

THEORY 
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A  method  of  finding  the  forward  solution  must  be  generated  first.  The  approach  taken  to 
the  defined  problem  assumes  that  a  ray  tracing  is  applicable  and,  as  previously  mentioned, 
that  the  impact  signals  move  out  in  all  half  space  directions.  It  is  instructive  to  look  at  a 
series  of  problems  of  increasing  difficulty  to  understand  the  complications  added  to  the 


F  i  gure  3.  Ti  me  versus  radial  di  stance  curve,  Fi  gure  4.  Geometry  for  Sen  sor  Pair 

isotropic  location  problem.  These  are  the  thick  isotropic  plate  problem,  the  multilayer 
thick  isotropic  plate  problem,  the  thick  anisotropic  plate  problem,  and  the  multilayer  thick 
anisotropic  plate  problem.The  thick  isotropic  plate  problem  is  directly  solvable 
analytically.  An  alternate  method,  however,  may  be  used.  It  is  noticed  that,  once  the 
thickness  is  fixed,  then  the  propagation  time  is  a  function  of: 

r  =  ^(jc,  -xj  +  (y,  -y,): 

only.  This  information  is  used  to  generate  a  curve  of  the  propagation  time,  tr  versus  the 
radial  distance,  r,  from  the  impact  axis  to  the  receipt  location.  This  t(r)  curve,  figure  3,  has  a 
non-zero  value  of  t  for  a  zero  radial  propagation  distance.  Also,  the  slope  of  the  curve  is  not 
constant  even  though  the  material  is  isotropic. 

Once  the  t(r)  function  is  found,  then  the  inverse  function  r(t)  can  be  generated.  If  the 
transit  time  of  the  impact  signal  to  two  sensors  were  known  then  the  backward  problem 
could  be  solved  by  using  r(t)  and  using  triangulation.  Additional  information  is  required 
since  the  initial  impact  time  is  unknown.  What  may  be  found  from  the  relative  arrival  times 
for  the  two  sensors  is  the  locus  of  possible  impact  points. 

The  constant  time  difference  curves  are  constructed  in  the  sensor  pair  local  coordinate 
system  shown  in  figure  4  in  the  parametric  form: 

*=jr'(s),  y  =s 

The  curve  for  a  particular  time  difference  is  created  by  assuming  a  propagation  time  to 
one  of  the  sensors.  The  minimal  value  of  this  propagation  time,  t(ri),  is  determined  by  the 


distances  2c*  between  sensors.  The  point  of  possible  impact  on  the  line  connecting  the 
sensors  is: 

where  r/  is  the  minimum  value  of  r/. 


Figure  5.  Family  of  constant  lime  difference  curves  Figure  6.  Relationship  of  Sensor  Local 


for  a  Sensor  Pair. 


Coordinates  to  Global  Coordinates. 


Once  r*  is  found,  other  points  on  the  curve  may  be  found  by  gradually  increasing  r/  and 
solving: 

r7  =HH-4)) 


When  the  curves  are  transformed  into  the  global  coordinate  system  (figure  6)  the 
parametric  equations  become: 
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By  adding  a  third  sensor  there  are  three  unique  sensor  pairs.  Two  of  the  sensor  pairs  may 
be  used  for  the  solution  and  the  intersection  of  their  possible  impact  curves  is  the  impact 
point.  If  timing  information  is  without  error  then  the  third  sensor  pair  is  redundant 
information  because  all  of  the  curves  will  intersect  at  the  same  point.  Otherwise,  the  third 
sensor  may  be  used  to  improve  the  estimate  of  the  impact  point.  If  a  point  is  picked  on  each 
curve  as  an  estimate  of  the  actual  impact  point  one  measure  of  the  merit  of  the  estimate 
would  be  the  sum  of  the  distances  between  the  estimate  points  on  these  curves.  This  leads  to 
an  objective  function; 

s  =  £  S  ((*<■  fo )  -  ■ x;  (*  + (ft  (*ihyj  ))* ) 

e=L  _^r+l 


to  be  minimized.  The  impact  point  estimates  on  the  curves  which  minimize  this  would  be 
the  point  on  each  curve  which  was  the  best  impact  point  estimate.  A  global  best  estimate 
can  be  found  from  taking  the  centroid  of  the  set  of  curve  estimate  points.  This  approach 
allows  additional  sensors  to  be  added  for  redundancy  in  a  natural  manner.  For  N  sensors 
there  will  be  N(N-l)/2  unique  sensor  pairs.  Using  an  explicitly  differentiable 
representation  for  these  curves  allows  the  use  of  a  gradient  search  method  to  minimize  the 
objective  function.. 

Adding  a  second  isotropic  layer  to  the  problem  requires  the  calculation  of  an  intercept 
point  at  the  interface  between  the  layers.  The  transmission  angle  of  the  mode  must  also  be 
found.  A  t(r)  curve  may  be  calculated  for  multiple  isotropic  layers  by  assuming  an  initial 
propagation  angle.  Using  this  propagation  angle,  the  radial  propagation  distance  and  transit 
time  are  calculated  for  the  current  layer  and  the  transmission  angle  into  the  next  layer  is 
calculated  using  Snell’s,  The  transmission  angle  is  used  for  calculation  the  radial 
propagation  distance  for  the  next  layer  and  the  transmission  angle  is  found  for  the 
subsequent  angle.  This  process  is  repeated  until  the  final  layer  is  processed.  A  running  total 
of  the  radial  distance  and  propagation  time  is  kept.  By  using  a  range  of  initial  angles,  this 
‘shooting’  method  allows  the  construction  of  a  series  of  points  on  the  t(r)  curve.  If  at  any 
stage  reflection  occurs  at  the  interface  then  the  bottom  of  the  assembly  cannot  be  reached. 
The  largest  previously  reached  radial  distance  is  considered  to  be  the  ‘radius  of 
observability’.  The  maximum  ‘shooting’  angle  used  is  capped  at  80  deg.  for  practical 
reasons,  even  though  90  deg  is  theoretically  possible. 

Consideration  of  the  anisotropic  layer  requires  the  introduction  of  the  concepts  of  the 
ChristofFel  equations,  phase  slowness,  phase  velocity,  and  group  velocity.  In  an  anisotropic 
material  the  propagation  velocity  of  a  mode  is  dependent  upon  the  propagation  direction. 
There  are  three  independent  modes,  as  with  an  isotropic  medium,  but  they  are  no  longer 
pure  modes.  This  means  that,  except  for  certain  special  directions,  the  modes  are  not  either 
parallel  or  perpendicular  to  the  direction  of  maximal  phase  velocity.  The  direction  of  energy 
flux,  the  group  velocity  direction,  is  not  parallel  to  phase  velocity  as  with  an  isotropic 
medium.  The  phase  slowness  is  the  reciprocal  of  the  phase  velocity. 

The  Christoffel  equations  [!  ]  allow  the  calculation  of  the  phase  slowness  for  a  particular 
set  of  direction  cosines.  This  eigenvalue  problem  yields  three  mode  phase  slowness 
magnitudes(eigen values)  and  directions{eigenvectors).  Evaluation  over  all  propagation 
directions  gives  a  set  of  three  closed  surfaces  as  in  figure  7.  One  surface  normally  is  most 
closely  aligned  with  the  propagation  direction,  has  the  highest  propagation  velocity,  and  is 
known  as  the  quasilongitudinal  mode.  There  are  two  other  modes  known  as  quasishear.  The 
phase  slowness  surfaces  are  spherical  in  an  isotropic  material.  In  a  transversely  isotropic 
material  they  are  axisymmetric  about  the  symmetry  axis.  Thus  only  the  curves  for  any  cross 
section  passing  through  the  symmetry  axis  need  be  calculated. 

The  group  velocity  direction[2]  is  normal  to  the  phase  slowness  surface  and  its 
projection  on  the  phase  velocity'  direction  vector  has  a  magnitude  equal  to  the  phase  velocity 
magnitude  as  shown  in  figure  8.  The  perceived  propagation  direction  in  an  anisotropic 
medium  is  the  group  velocity.  The  phase  velocity  and  group  velocity  surfaces  may  be 
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constructed  once  the  phase  slowness  surface  is  obtained.  The  t(r)  curve  for  the  layer  may  be 
found  using  the  ‘shooting  method’  with  the  group  velocity. 


Finally,  the  ease  of  multiple  anisotropic  layers  requires  the  use  of  the  generalized  Snell’s 
law  as  shown  in  figure  9.  This  results  from  the  requirement  that  the  transverse  component  of 
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for  transversely  isotropic  material. 


Figure  S.  Relationship  between  Phase 
Slowness  and  Group  Velocity  Directions . 


phase  slowness  be  equal  at  an  interface  between  two  dissimilar  materials.  The  procedure 
for  ray  tracing  through  the  assembly  is  to  assume  an  initial  propagation  angle  for  the 
group  velocity  through  the  first  layer.  The  radial  propagation  and  transit  time  through  the 
layer  are  then  calculated.  The  point  on  the  phase  slowness  surface  corresponding  to  the 
group  velocity  direction  is  located.  The  generalized  Snell’s  law  is  used  to  locate  the 
corresponding  point  on  the  phase  slowness  surface  of  the  next  layer.  The  group  velocity 
in  the  next  layer  is  found  from  its*  association  with  the  point  on  the  phase  slowness 
surface.  This  process  is  repeated  until  the  bottom  of  the  last  layer  is  reached. 


IMPLEMENTATION 


This  procedure  has  been  implemented  as  an  application  in  MATLAB®.  The  phase 
slowness  curves  of  the  materials  in  the  model  are  generated  by  solving  the  Christoffel 
Equations  in  range  of  propagation  directions  in  the  XZ  plane.  The  x  and  z  components  of  the 
points  on  the  curves  are  then  parameterized  in  terms  of  phi  using  nonlinear  regression.  This 
allows  the  tangent  and  normal  vectors  to  be  found  from  the  derivatives  of  the  parameterized 
coordinate  equations. 

Once  the  regression  curves  for  the  phase  slowness  components  are  generated  the 
phase  velocity  curves  can  be  calculated  from  the  inverse  of  the  phase  slowness 
magnitude.  The  group  velocity  curves  are  then  found  by  use  of  the  normal  to  the  phase 
slowness  surface  and  the  phase  velocity  magnitude.  This  information  is  stored  in  the  form 
of  group  velocity  magnitude  and  direction. 
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Figure  9,  Generalized  Snell's  Law 


Figure  10.  Impact  Location  Algorithm  output  for  Figure  1 1.  Close  Up  of  impact  Location 
Observable  Region.  RegionforAl  gerilhm  output 

The  group  velocity  information  for  each  layer  is  used  to  determine  the  radial 
propagation  distance  and  transit  time.  Linear  interpolation  is  used  when  the  propagation 
angle  falls  between  known  values.  If  there  is  more  than  one  layer  the  group  velocity  and 
phase  slowness  information  of  touching  layers  is  used  together  for  solution  of  the 
generalized  Snell’s  law  relation.  Once  again  linear  interpolation  is  used. 

Upon  finishing  the  time  versus  radial  distance  curve,  the  constant  time  of  flight 
difference  curves  are  generated  for  each  sensor  pair.  This  is  performed  by  observing 
which  sensor  has  the  smallest  absolute  time  of  arrival.  This  indicates  which  sensor  is 
closest  to  the  intercept  of  the  constant  time  of  flight  curve  and  the  sensor  centerline,  Ibis 
intercept  point  is  found  using  bisection.  The  other  points  on  the  constant  time  of  flight 
curve  within  the  radius  of  observability  are  then  calculated  by  using  linear  interpolation  to 
solve  the  relation  between  rj  and  rj  given  previously. 

These  curves  are  approximated  over  the  radius  of  observability  using  Chebyshev 
polynomials [3].  The  Chebyshev  polynomials  are  then  converted  to  regular  polynomials 
using  the  method  outlined  in  [4],  The  results  of  the  approximation  are  shown  in  figure  5. 
The  regular  polynomials  are  then  converted  from  the  local  sensor  pair  coordinates  into 
global  coordinates  using  the  transformation  previously  given. 

These  regular  polynomial  curves  are  used  to  construct  the  objective  function.  The  set  of 
coordinate  parameter  values  that  minimize  the  objective  function  are  found  using  the 
Newton-Kantorovich  algorithm^  ]  as  follows.  An  initial  guess  of  the  coordinate  parameter 
for  each  curve  is  made.  The  Hessian  matrix  and  gradient  of  the  objective  function  is  then 
evaluated  at  the  current  coordinate  parameter  values.  The  Hessian  matrix  inverse  and  the 
gradient  vector  are  then  used  to  create  an  updated  coordinate  parameter  estimate.  This  is 
continued  until  convergence  criteria  are  met  or  divergence  is  detected.  If  convergence 
occurs,  then  the  final  estimate  of  the  impact  location  is  found  by  evaluation  the  coordinates 
of  each  sensor  pair  curve  from,  the  current  coordinate  parameter  values  and  then  finding  die 
centroid  of  this  set  of  values  as  illustrated  in  figure  1 1 . 
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RESULTS 


The  method  was  evaluated  using  simulated  timing  data  for  layers  of  isotropic  and 
anisotropic  materials.  It  was  also  tested  against  experimental  data  gathered  from  impacts 
on  an  aluminum  block.  Figures  10  and  1 1  show  graphical  output  from  the  program  for  a 
typical  lost  data  ease. 

These  figures  show  the  sensor  locations,  the  actual  impact  point,  the  constant  time 
difference  curves  generated  for  the  sensor  pairs  from  the  timing  data,  the  evaluation  of  the 
initial  estimates  on  the  curves,  the  intermediate  estimates,  the  converged  impact  location 
estimates,  and  the  final  impact  location  estimate.  Figure  10  shows  these  entities  over  the 
entire  observable  area.  Figure  1 1  shows  a  close-up  of  the  impact  location  area, 

CONCLUSIONS 

A  method  of  finding  impact  location  by  using  timing  data  from  sensors  residing  on 
the  backside  of  an  assembly  composed  of  multiple  layers  of  transversely  isotropic 
material  has  been  formulated.  This  method  works  by  solving  a  Model  based  Inverse 
problem  where  the  propagation  time  through  a  multilayer  transversely  isotropic  layer 
assembly  is  the  forward  problem.  The  method  has  been  implemented  in  software  and 
tested  against  simulated  data  for  the  multilayer  transversely  isotropic  case  and  actual  data 
for  the  single  layer  isotropic  case, 
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